Demo of the L-Wigner Distribution and the Polynomial WV from STFT & SM

last update: 2/15 (2018)

Based on algorithms in Boashash+ 2015 p349 (6.2.17) & p350 (6.2.19)



In [1]:

    
using PyPlot
import DSP



In [2]:

    
include("../juwvid.jl")









    Out[2]:





juwvid

multi frequency-components data (Boashash+15,p346,Example 6.2.3)



In [3]:

    
nsample=512
t,x=sampledata.genmultifm623(nsample);
PyPlot.plot(t,x)
PyPlot.ylabel("y")
PyPlot.xlabel("t")









    












    Out[3]:





PyObject Text(0.5,24,'t')

S-method (constant Lp)



In [4]:

    
sm=smethod.tfrsm(x,NaN,4,NaN,2);









    



Single S-method
Use fft.



In [5]:

    
a=juwplot.wtfrshow(sm,t[2]-t[1],t[1],t[end],NaN,NaN,1.0)
PyPlot.xlabel("time")
PyPlot.ylabel("frequency")
PyPlot.ylim(0,65)









    












    Out[5]:





(0, 65)

L-Wigner Distribution (L=2)



In [6]:

    
#alias free sm
afwv=smethod.tfrsm(x,NaN,4,NaN,2)
#L wigner distribution (L=2)
tfrlw2=lwigner.tfrlw2L(afwv,4);
a=juwplot.wtfrshow(tfrlw2,t[2]-t[1],t[1],t[end],NaN,NaN,1.0)
PyPlot.xlabel("time")
PyPlot.ylabel("frequency")
PyPlot.ylim(0,65)









    












    



Single S-method
Use fft.






    Out[6]:





(0, 65)

polynomial Wigner Ville distribution



In [7]:

    
trfpo=polywv.tfrpowv(x,NaN,NaN,NaN,2,8);
a=juwplot.wtfrshow(trfpo,t[2]-t[1],t[1],t[end],NaN,NaN,1.0)
PyPlot.xlabel("time")
PyPlot.ylabel("frequency")
PyPlot.ylim(0,65)









    












    



Single S-method
Use fft.






    Out[7]:





(0, 65)

comparison



In [8]:

    
tfrst=stft.tfrstft(x);









    



Use fft.



In [9]:

    
fig=PyPlot.figure()
ax = fig[:add_subplot](2,2,1)
a=juwplot.wtfrshow(abs.(tfrst).^2,t[2]-t[1],t[1],t[end],NaN,NaN,0.7*2)
PyPlot.ylabel("frequency")
PyPlot.title("Spectrogram")
PyPlot.ylim(0,65)
ax = fig[:add_subplot](2,2,2)
a=juwplot.wtfrshow(abs.(sm),t[2]-t[1],t[1],t[end],NaN,NaN,0.7*2)
PyPlot.title("S-method")
PyPlot.ylim(0,65)
ax = fig[:add_subplot](2,2,3)
a=juwplot.wtfrshow(abs.(tfrlw2),t[2]-t[1],t[1],t[end],NaN,NaN,0.7*2)
PyPlot.xlabel("time")
PyPlot.ylabel("frequency")
PyPlot.title("L-wigner")
PyPlot.ylim(0,65)
ax = fig[:add_subplot](2,2,4)
a=juwplot.wtfrshow(abs.(trfpo),t[2]-t[1],t[1],t[end],NaN,NaN,0.7*2)
PyPlot.xlabel("time")
PyPlot.title("Polynomial WV")
PyPlot.ylim(0,65)









    












    Out[9]:





(0, 65)

Comparison (noisy case)



In [10]:

    
using Distributions
d = Normal()
xnoise=x+rand(d,nsample)*std(x)*0.5
PyPlot.plot(t,xnoise)









    












    Out[10]:





1-element Array{PyCall.PyObject,1}:
 PyObject <matplotlib.lines.Line2D object at 0x7feabfae20f0>



In [11]:

    
tfrst=stft.tfrstft(xnoise);
#alias free sm
sm=smethod.tfrsm(xnoise,NaN,4,NaN,2)
#L wigner distribution (L=2)
tfrlw2=lwigner.tfrlw2L(sm,4);
trfpo=polywv.tfrpowv(xnoise,NaN,NaN,NaN,2,8);









    



Use fft.
Single S-method
Use fft.
Single S-method
Use fft.



In [12]:

    
fig=PyPlot.figure()
ax = fig[:add_subplot](2,2,1)
a=juwplot.wtfrshow(abs.(tfrst).^2,t[2]-t[1],t[1],t[end],NaN,NaN,0.7*2)
PyPlot.ylabel("frequency")
PyPlot.title("Spectrogram")
PyPlot.ylim(0,65)
ax = fig[:add_subplot](2,2,2)
a=juwplot.wtfrshow(abs.(sm),t[2]-t[1],t[1],t[end],NaN,NaN,0.7*2)
PyPlot.title("S-method")
PyPlot.ylim(0,65)
ax = fig[:add_subplot](2,2,3)
a=juwplot.wtfrshow(abs.(tfrlw2),t[2]-t[1],t[1],t[end],NaN,NaN,0.7*2)
PyPlot.xlabel("time")
PyPlot.ylabel("frequency")
PyPlot.title("L-wigner")
PyPlot.ylim(0,65)
ax = fig[:add_subplot](2,2,4)
a=juwplot.wtfrshow(abs.(trfpo),t[2]-t[1],t[1],t[end],NaN,NaN,0.7*2)
PyPlot.xlabel("time")
PyPlot.title("Polynomial WV")
PyPlot.ylim(0,65)









    












    Out[12]:





(0, 65)